

* Title: Behind political affiliation: how moral values, identity politics, and party loyalty have affected COVID-19 vaccination

*************************************************************************************************
***    PRELIMINARY SETTING   ***
*************************************************************************************************
*clear any data currently in Stata's memory 
 clear all

*set command interpreter to version 17
version 17

*tell Stata to run the entire do-file without pausing
set more off


 
***********************************************************************************
**# 1° RC: Baseline with FIPS-FE
***********************************************************************************


qui{
. use ".\anti-vaxxer dataset - temporary.dta", replace

xtset num_fips date
 global predictors Unemployment_rate_2021 Median_HHIncome2020  poverty_rate prop_lessthan_HSedu pop_density internet_subsc ep_age65 ep_minrty health_insurance_rate staffedallbedsper1000people

 
 
**# S3 
// Table Political partisanship and COVID-19 vaccination rates with county fixed effects.

* Regression 1. 
  eststo clear
 eststo baseline1: reghdfe series_complete_pop_pct i.rep_majority_2020##i.qdate L7.(total_cases_rate total_deaths_rate), a(mdate fips i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 


* Regression 2. 
 eststo baseline2: reghdfe series_complete_pop_pct i.partisanship_2020##i.qdate L7.(total_cases_rate total_deaths_rate), a(mdate fips i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 


* Regression 3.
 eststo baseline3: reghdfe series_complete_pop_pct i.rooted_partisanship##i.qdate L7.(total_cases_rate total_deaths_rate), a(mdate fips i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 


* Regression 4.
 eststo baseline4: reghdfe series_complete_pop_pct i.trumpist_county##i.qdate L7.(total_cases_rate total_deaths_rate)  if rep_majority_2020==1, a(mdate fips i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 

}
 
 


 esttab baseline1 baseline2 baseline3 baseline4 using ".\S3 table.tex", stats(N r2,  labels("Observations" "R-squared" )) compress replace b(3) se(3) nonotes nogaps  starlevels(* 0.10 ** 0.05 *** 0.01) mtitle("Model (1)" "Model (2)" "Model (3)" "Model (4)") drop( 2**.qdate ***244.qdate 0.rep_majority_2020*** 1.partisanship_2020*** *.partisanship_2020 1.rooted_partisanship*** *.rooted_partisanship 0.trumpist_county*** 1.trumpist_county) title("Political partisanship and Covid-19 vaccination rates with county fixed effects") longtable ///
 addnotes("The dependent variable is Fully vaccinated people (%). Models 1-4 include FIPS, month and state*month fixed effects. Model 4 is conducted only on counties where Republicans obtained the majority of votes in the 2020 presidential election. Q1, Q2, Q3, and Q4 refer to the periods: Jan 1–Mar 31, Apr 1–Jun 30, Jul 1–Sep 30, and Oct 1–Dec 31, respectively. All models use two-way robust standard errors clustered by state and week. Robust standard errors in parentheses. *** p<0.01, ** p<0.05, * p<0.1.") label
 
 
 
 
 ***********************************************************************************
**# 2° RC:  Analysis with time-centered data
***********************************************************************************

 . use ".\anti-vaxxer dataset - temporary.dta", replace

 mdesc series_complete_yes series_complete_pop_pct 
 
by fips (date), sort: gen campaign_started = sum(series_complete_pop_pct) > 0


by fips (date), sort: gen days_centered = sum(campaign_started[_n]) if campaign_started ==1


order campaign_started days_centered, after(qdate)

drop if campaign_started == 0
sum days_centered

* Gen log dep. variables
  gen ln_vaccinations = ln(series_complete_yes)
 label var ln_vaccinations "ln(Total number of fully vaccinated people)"
 

**# S4 Table 
// Political partisanship and COVID-19 vaccination rates with time-centered vaccination data.

 xtset num_fips days_centered
 
global predictors Unemployment_rate_2021 Median_HHIncome2020  poverty_rate prop_lessthan_HSedu pop_density internet_subsc ep_age65 ep_minrty health_insurance_rate staffedallbedsper1000people



* Regression 1. 
  eststo clear
 eststo baseline1: reghdfe series_complete_pop_pct i.rep_majority_2020##i.qdate L7.(total_cases_rate total_deaths_rate) $predictors, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 
estadd local controlsvar "Yes"


* Regression 2. 
 eststo baseline2: reghdfe series_complete_pop_pct i.partisanship_2020##i.qdate L7.(total_cases_rate total_deaths_rate) $predictors, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 
estadd local controlsvar "Yes"

 
* Regression 3.
 eststo baseline3: reghdfe series_complete_pop_pct i.rooted_partisanship##i.qdate L7.(total_cases_rate total_deaths_rate) $predictors, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 
estadd local controlsvar "Yes"


* Regression 4.
 eststo baseline4: reghdfe series_complete_pop_pct i.trumpist_county##i.qdate L7.(total_cases_rate total_deaths_rate) $predictors if rep_majority_2020==1, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 
estadd local controlsvar "Yes"


 
esttab baseline1 baseline2 baseline3 baseline4 using ".\S4 table.tex", stats(N r2 controlsvar,  labels("Observations" "R-squared" "Health and economic controls")) compress replace b(3) se(3) nonotes nogaps  starlevels(* 0.10 ** 0.05 *** 0.01) mtitle("Model (1)" "Model (2)" "Model (3)" "Model (4)") drop($predictors 2**.qdate ***244.qdate 0.rep_majority_2020*** 1.partisanship_2020*** 1.rooted_partisanship*** 0.trumpist_county***) ///
 title("Political partisanship and Covid-19 vaccination rates with time-centered vaccination campaign") longtable ///
 addnotes("The dependent variable is Fully vaccinated people (%). Models 1-4 include state, month and state*month fixed effects. Model 4 is conducted only on counties where Republicans obtained the majority of votes in the 2020 presidential election. Health and economic controls include: unemployment rate, poverty rate, median income, proportion of people with less than HS diploma, population density, proportion of people with internet subscription, proportion of people older than 65 years old, proportion of minority, proportion of people with health insurance, health insurance rate number of hospital beds of all types per 1000 people. Q1, Q2, Q3, and Q4 refer to the periods: Jan 1–Mar 31, Apr 1–Jun 30, Jul 1–Sep 30, and Oct 1–Dec 31, respectively. All models use two-way robust standard errors clustered by state and week. Robust standard errors in parentheses. *** p<0.01, ** p<0.05, * p<0.1.") label
 
  

  

 ***************************************************************************************
**# 3° RC: Political partisanship variables in 2016
*************************************************************************************** 
 use ".\anti-vaxxer dataset - temporary.dta", replace

xtset num_fips date
 global predictors Unemployment_rate_2021 Median_HHIncome2020  poverty_rate prop_lessthan_HSedu pop_density internet_subsc ep_age65 ep_minrty health_insurance_rate staffedallbedsper1000people


 label var rep_majority_2016 "Republican county (2016)"
 
**# S5 Table
// Political partisanship (in 2016) and Covid-19 vaccination rates
 
* Regression 1. 
  eststo clear
 eststo baseline1: reghdfe series_complete_pop_pct i.rep_majority_2016##i.qdate L7.(total_cases_rate total_deaths_rate) $predictors, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 
estadd local controlsvar "Yes"


* Regression 2. 
 eststo baseline2: reghdfe series_complete_pop_pct i.partisanship_2016##i.qdate L7.(total_cases_rate total_deaths_rate) $predictors, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 
estadd local controlsvar "Yes"

 
* Regression 3.
 eststo baseline3: reghdfe series_complete_pop_pct i.rooted_partisanship2016##i.qdate L7.(total_cases_rate total_deaths_rate) $predictors, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 
estadd local controlsvar "Yes"


* Regression 4.
 eststo baseline4: reghdfe series_complete_pop_pct i.trumpist_county2016##i.qdate L7.(total_cases_rate total_deaths_rate) $predictors if rep_majority_2020==1, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 
estadd local controlsvar "Yes"



 
esttab baseline1 baseline2 baseline3 baseline4 using ".\S5 table.tex", stats(N r2 controlsvar,  labels("Observations" "R-squared" "Health and economic controls")) compress replace b(3) se(3) nonotes nogaps  starlevels(* 0.10 ** 0.05 *** 0.01) mtitle("Model (1)" "Model (2)" "Model (3)" "Model (4)") drop($predictors 2**.qdate ***244.qdate 0.rep_majority_2016*** 1.partisanship_2016*** 1.rooted_partisanship2016*** 0.trumpist_county2016***) ///
 title("Political partisanship and Covid-19 vaccination rates (in 2016)") longtable ///
 addnotes("The dependent variable is Fully vaccinated people (%). Models 1-4 include state, month and state*month fixed effects. Model 4 is conducted only on counties where Republicans obtained the majority of votes in the 2020 presidential election. Health and economic controls include: unemployment rate, poverty rate, median income, proportion of people with less than HS diploma, population density, proportion of people with internet subscription, proportion of people older than 65 years old, proportion of minority, proportion of people with health insurance, health insurance rate number of hospital beds of all types per 1000 people. Q1, Q2, Q3, and Q4 refer to the periods: Jan 1–Mar 31, Apr 1–Jun 30, Jul 1–Sep 30, and Oct 1–Dec 31, respectively. All models use two-way robust standard errors clustered by state and week. Robust standard errors in parentheses. *** p<0.01, ** p<0.05, * p<0.1.") label
  
  
  
  
  
  ***************************************************************************************
**# 4° RC: Testing for multicollinearity
*************************************************************************************** 
 use ".\anti-vaxxer dataset - temporary.dta", replace

xtset num_fips date
 global predictors Unemployment_rate_2021 Median_HHIncome2020  poverty_rate prop_lessthan_HSedu pop_density internet_subsc ep_age65 ep_minrty health_insurance_rate staffedallbedsper1000people

// Generate a new variable marking the last day of each month
gen month = mofd(date)
format month %tm // OPTIONAL IF YOU DON'T NEED THIS VARIABLE FOR OTHER PURPOSES
by fips month (date), sort: keep if _n == _N

xtset num_fips mdate
 
 
 
**# S6 Table
// VIF for political and moral variables.

// Define the variables
global political_vars rep_majority_2020 partisanship_2020 rooted_partisanship trumpist_county 
global moral_vars relative_comm_1518 values_communal1518_shrink 

	
// Loop through political and moral variables
foreach pol_var in $political_vars {
		foreach moral_var in $moral_vars {

			// Run the regression
			qui: regress series_complete_pop_pct `pol_var' `moral_var' ///
				total_cases_rate total_deaths_rate $predictors ///
				i.qdate##i.state_ID, cluster(state_ID)

			// Add the VIF values to the current model estimation
			qui: estat vif
			qui: return list
			
		  // Save and round VIF values
        forvalues i = 1/2 {
            // Round the VIF value to 2 decimal places
            local vif_value = round(`r(vif_`i')', 0.01)
            di "Variable: `r(name_`i')' VIF: `vif_value'"
			}
			di "                            "
		}
}




***************************************************************************************
**# 5° RC: Testing for FDR
*************************************************************************************** 
// S7 Table. Sharpened q-values for primary regression model specifications.

do "code3_fdr_sharpened_qvalues.do"

 save ".\qvals_all.dta", replace
 
 
  use ".\anti-vaxxer dataset - temporary.dta", replace

xtset num_fips date
 global predictors Unemployment_rate_2021 Median_HHIncome2020  poverty_rate prop_lessthan_HSedu pop_density internet_subsc ep_age65 ep_minrty health_insurance_rate staffedallbedsper1000people

 gen year = 2021
 replace year =2022 if qdate ==tq(2022q1)
  replace year =2022 if qdate ==tq(2022q2)
 

 
***************************************************************************************
**# 6° RC: Testing for commuting zone and quarters FE
*************************************************************************************** 
// To account for changing in information supply and consumption across communities
// Data for commuting zones can be retrieved from USDA: https://www.ers.usda.gov/data-products/commuting-zones-and-labor-market-areas 


 use ".\anti-vaxxer dataset - temporary.dta", replace



  merge m:1 fips using ".\cz00_eqv_v1.dta"
  
  drop if _merge ==2 
  label var CommutingZoneID2000 "Commuting zone"
 
 xtset num_fips date
 
 **# S8 Table

* Regression 1. 
  eststo clear
 eststo baseline1: reghdfe series_complete_pop_pct i.rep_majority_2020##i.qdate $predictors L7.(total_cases_rate total_deaths_rate), a(mdate state i.CommutingZoneID2000#i.qdate) cluster(i.state_ID i.wdate) 


* Regression 2. 
 eststo baseline2: reghdfe series_complete_pop_pct i.partisanship_2020##i.qdate $predictors L7.(total_cases_rate total_deaths_rate), a(mdate state i.CommutingZoneID2000#i.qdate) cluster(i.state_ID i.wdate) 


* Regression 3.
 eststo baseline3: reghdfe series_complete_pop_pct i.rooted_partisanship##i.qdate $predictors L7.(total_cases_rate total_deaths_rate), a(mdate state i.CommutingZoneID2000#i.qdate) cluster(i.state_ID i.wdate) 


* Regression 4.
 eststo baseline4: reghdfe series_complete_pop_pct i.trumpist_county##i.qdate $predictors L7.(total_cases_rate total_deaths_rate)  if rep_majority_2020==1, a(mdate state i.CommutingZoneID2000#i.qdate) cluster(i.state_ID i.wdate) 
 
 
 
 esttab baseline1 baseline2 baseline3 baseline4 using ".\S8 table.tex", stats(N r2,  labels("Observations" "R-squared" )) compress replace b(3) se(3) nonotes nogaps  starlevels(* 0.10 ** 0.05 *** 0.01) mtitle("Model (1)" "Model (2)" "Model (3)" "Model (4)") drop( 2**.qdate ***244.qdate 0.rep_majority_2020*** 1.partisanship_2020*** *.partisanship_2020 1.rooted_partisanship*** *.rooted_partisanship 0.trumpist_county*** 1.trumpist_county) title("Political partisanship and Covid-19 vaccination rates with CZ fixed effects") longtable ///
 addnotes("The dependent variable is Fully vaccinated people (%). Models 1-4 include state, month and commuting zone*quarter fixed effects. Model 4 is conducted only on counties where Republicans obtained the majority of votes in the 2020 presidential election. Q1, Q2, Q3, and Q4 refer to the periods: Jan 1–Mar 31, Apr 1–Jun 30, Jul 1–Sep 30, and Oct 1–Dec 31, respectively. All models use two-way robust standard errors clustered by state and week. Robust standard errors in parentheses. *** p<0.01, ** p<0.05, * p<0.1.") label
 

 
 
 
 ***************************************************************************************
**# 7° RC: Testing for the effect of fairness perception of vaccination mandates
*************************************************************************************** 
// Robustness check not presented in the final draft
// Relevant paper: https://jamanetwork.com/journals/jamanetworkopen/fullarticle/2774317?utm_source=chatgpt.com
 
 
use ".\anti-vaxxer dataset - temporary.dta", replace

 
xtset num_fips date

gen copartisanship_dem =0
replace copartisanship_dem = 1 if rooted_partisanship2000_20==1 & governor_party_num==0
 
 
 
// Generate indicator variable for states with vaccination exemptions
// List of states provided in: https://jamanetwork.com/journals/jama-health-forum/fullarticle/2797733?utm_source=chatgpt.com

* Generate a variable for mandate group states
gen mandate_group = 0
replace mandate_group = 1 if inlist(State, "CA", "CO", "CT", "DC", "ME") | inlist(State, "MA", "NV", "NM", "NY") | inlist(State, "NC", "OR", "RI", "WA")

* Generate a variable for comparison group states
gen comparison_group = 0
replace comparison_group = 1 if inlist(State, "DE", "HI", "IL", "KY", "LA") | inlist(State, "MD","MN", "NJ", "OH", "PA") |inlist(State, "VA", "VT", "WV", "WI")

* Generate a variable for states with a prohibition of mandates (remaining states)
gen prohibition_group = 0
replace prohibition_group = 1 if !(mandate_group | comparison_group)


*** Baseline analysis on states with a prohibition of mandates

 * Regression 1. 
  eststo clear
 eststo baseline1: reghdfe series_complete_pop_pct i.rep_majority_2020##i.qdate L7.(total_cases_rate total_deaths_rate) if prohibition_group==1, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 


* Regression 2. 
 eststo baseline2: reghdfe series_complete_pop_pct i.partisanship_2020##i.qdate L7.(total_cases_rate total_deaths_rate) if prohibition_group==1, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 


* Regression 3.
 eststo baseline3: reghdfe series_complete_pop_pct i.rooted_partisanship##i.qdate L7.(total_cases_rate total_deaths_rate) if prohibition_group==1, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 
 
 * Regression 4.
  eststo baseline4: reghdfe series_complete_pop_pct i.trumpist_county##i.qdate $predictors L7.(total_cases_rate total_deaths_rate)  if prohibition_group==1, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 

 
 
 esttab baseline1 baseline2 baseline3 baseline4 using ".\S_extra table.rtf", stats(N r2,  labels("Observations" "R-squared" )) compress replace b(3) se(3) nonotes nogaps  starlevels(* 0.10 ** 0.05 *** 0.01) mtitle("Model (1)" "Model (2)" "Model (3)" "Model (4)") drop( 2**.qdate ***244.qdate 0.rep_majority_2020*** 1.partisanship_2020*** *.partisanship_2020 1.rooted_partisanship*** *.rooted_partisanship) title("Political partisanship and Covid-19 vaccination rates (states that issued mandate prohibitions)") longtable ///
 addnotes("The dependent variable is Fully vaccinated people (%). Models 1-4 include state, month and state*month fixed effects. Q1, Q2, Q3, and Q4 refer to the periods: Jan 1–Mar 31, Apr 1–Jun 30, Jul 1–Sep 30, and Oct 1–Dec 31, respectively. Model 4 is conducted only on counties where Republicans obtained the majority of votes in the 2020 presidential election. All models use two-way robust standard errors clustered by state and week. Robust standard errors in parentheses. *** p<0.01, ** p<0.05, * p<0.1.") label
 
 
 
 *** Baseline analysis on comparison group states
 
 * Regression 1. 
  reghdfe series_complete_pop_pct i.rep_majority_2020##i.qdate L7.(total_cases_rate total_deaths_rate) if comparison_group==1, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 


* Regression 2. 
reghdfe series_complete_pop_pct i.partisanship_2020##i.qdate L7.(total_cases_rate total_deaths_rate) if comparison_group==1, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 


* Regression 3.
  reghdfe series_complete_pop_pct i.rooted_partisanship##i.qdate L7.(total_cases_rate total_deaths_rate) if comparison_group==1, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 
 
 * Regression 4.
  reghdfe series_complete_pop_pct i.trumpist_county##i.qdate $predictors L7.(total_cases_rate total_deaths_rate)  if comparison_group==1, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 
 
 
 
 
 